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Abstract 

Consider a particle diffusing in a confined volume which is divided into two equal 
regions. In one region the diffusion coefficient is twice the value of the diffusion 
coefficient in the other region. Will the particle spend equal proportions of time in 
the two regions in the long term? Statistical mechanics would suggest yes, since 
the number of accessible states in each region is presumably the same. However, 
another line of reasoning suggests that the particle should spend less time in the 
region with faster diffusion, since it will exit that region more quickly. We demon- 
strate with a simple microscopic model system that both predictions are consistent 
with the information given. Thus, specifying the diffusion rate as a function of 
position is not enough to characterize the behaviour of a system, even assuming 
the absence of external forces. We propose an alternative framework for modelling 
diffusive dynamics in which both the diffusion rate and equilibrium probability 
density for the position of the particle are specified by the modeller. We introduce a 
numerical method for simulating dynamics in our framework that samples from the 
equilibrium probability density exactly and is suitable for discontinuous diffusion 
coefficients. 



1 Introduction 

Consider a particle diffusing in a two-dimensional box with reflecting boundary 
conditions. We show a portion of a simulated trajectory of such a system in Fig- 
ure[T| Suppose that in the left half of the box the particle diffuses with coefficient D\ 
and that in the right half of the box the particle diffuses with coefficient D2 = 2D\ . 
We assume that there are no external forces acting on the particle. Our question is: 
Does the particle spend an equal fraction of time on each side of the box in the long 
run? 

One answer is based on statistical mechanics. 

Statistical Mechanics Prediction: The particle will spend an equal propor- 
tion of time on each side of the box. 
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The justification for this prediction is the principle of statistical mechanics which 
states that "An isolated system in equilibrium is equally likely to be in any of its 
accessible states" (Reif 1965, p. 54). Since all states in the box are accessible, and 
there are an equal number of states on each side of the box, the particle should 
spend an equal proportion of its time on each side of the box. 

Another answer is based on the idea of rescaling time in one side of the box. 

Time-Change Prediction: The particle will spend less time on the side of the 
box where the diffusion coefficient is greater. 

The justification for this prediction is that, in the absence of any drift, faster dif- 
fusion is equivalent to time passing more quickly. This means that the periods of 
time the particle spends on the right side of the box will be shorter than those spent 
on the left side. Hence the total time the particle spends on the right hand side of 
the box will be less. We show how this prediction is a straightforward consequence 
of interpreting the particle's motion as drift-free diffusion, where we interpret the 
state-dependent diffusion coefficient using the Ito convention (Gardiner, 2004). We 
can write the equation for the particle's motion as 

dx = b(x)dB(t), (1.1) 

where x = (x\,X2) an d B is standard two-dimensional Brownian motion. (Equiva- 
lently we may write this equation as dx/dt = b{x)r\(t) where r\ is two-dimensional 
Gaussian white noise.) We specify b{x) = b\ for x\ < and b(x) = £2 for x\ > 0. 
Here b\ = \/2Di, i= 1,2 where Z), is the corresponding diffusion coefficient. We en- 
force reflecting boundary conditions at the four walls of the box. The (Ito-)Fokker- 
Planck equation for p(x,t), the probability density of the particle's position at time 
t,is (Gardiner 2004, p. 118) 

J^p(*,f) = • [V (b 2 (x)p(x,t))} = V • [V(Z>(*)p(*,f))]. 

The equilibrium density p eq (x) satisfies 

V (D(x)p eq (x)) = const. 

Reflecting boundary conditions for the diffusion correspond to Neumann boundary 
conditions (zero-flux) for the Fokker-Planck equation. With these boundary con- 
ditions, the unique equilibrium density is p eq (x) = C/D(x) for some constant C. 
Thus, since D2 = 2D\, the particle spends half as much time on the right side of the 
box as on the left. (We discuss the relation of our question to other interpretations 
of ( fTTTT ) in Section [4}) 

Neither the Statistical Mechanics Prediction nor the Time-Change Prediction are 
definitive. The Statistical Mechanics Prediction relies on the principle of equal 
probability of all accessible states, which needs to be independently justified for 
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Fig. 1. Simulation of particle diffusing within a box with reflecting boundary conditions. 
Points are generated by the Euler-Maruyama method X n+ \ = X n + y / 2hD(X n )N n , where 
h = 0.5, N„ are independent standard two-dimensional Gaussians, and D(x) is the state-de- 
pendent diffusion coefficient. Here D(x) = 1 on the left side of the box and D(x) =2 on the 
right side of the box. 

the mesoscopic level of description we are considering here. The Time-Change 
Prediction is not definitive since Ito stochastic differential equations are themselves 
mesoscopic models whose use can only be rigorously justified by showing how they 
arise as the coarse-scale limit of microscopic dynamics. Since the two predictions 
contradict each other, at least one of them must be wrong for any given physical 
system. The approach by which we resolve this apparent contradiction is to study a 
microscopic model of the box system we have described above and then see what 
proportion of the time the particle spends on each side in simulations of that sys- 
tem. If one prediction turned out to always be true for our model system, we could 
use our result as a baseline for investigating under which more general situations 
the prediction was still true. However, we will see that even for our simple system 
there are parameter choices that make the Statistical Mechanics Prediction correct 
and parameter choices that make the Time-Change Prediction correct. This demon- 
strates that there is no a priori reason to conclude that one prediction or another is 
correct, given only the diffusion coefficient D(x) for the system. 

In Section |2l we describe our model system, perform numerical experiments on it, 
and show that the result of the numerical experiments can be determined analyt- 
ically from the properties of the model system. The system we consider is a de- 
terministic Hamiltonian billiard system: the random Lorentz gas (Dettmann 2000). 
The system has two free parameters: disc radius and free volume fraction. There 
is a one-parameter family of values of these parameters that can generate the same 
effective diffusion coefficient. We will show how to choose these parameters to 
obtain arbitrarily good approximations to diffusive motion for the particle on each 
side of the box. Using the degree of freedom in the choice of parameters, we show 
that the equilibrium density of the particle is underdetermined by the diffusion co- 
efficient on each side of the region. For a given D\ and £>2, exploiting the flexibility 
in the parameter choice allows us to create systems in which either the Statistical 
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Mechanics Prediction or the Time-Change Prediction is correct. 

Our justification for studying a particular microscopic system is threefold. Firstly, 
there is a solid analytical understanding of the random Lorentz gas on which we can 
base our simulations. Secondly, our purpose is not to conclude that a particular style 
of mesoscopic modelling is always the correct one, but to show that postulating a 
state-dependent diffusion rate is not enough to fully specify mesoscopic behaviour. 
For this objective, it is enough to show that multiple mesoscopic behaviours are 
possible for a single class of simple models, as we do in this paper. Finally, we have 
resorted to a purely deterministic microscopic model, rather than a stochastic mi- 
croscopic model (such as Langevin dynamics or a random walk) to avoid concerns 
that the method used to introduce randomness at the microscopic level somehow 
biases the results at the mesoscopic level. This last point distinguishes our work 
from similar discussions of van Kampen (2007, Sec. XI.3) and Othmer and Stevens 
(1997, Sec. 2). These authors show that spatially inhomogeneous random walks, 
with natural choices of parameters, can lead to either prediction holding true at the 
mesoscopic level. 

Similarly, the conclusions of our Section [2] closely parallel those of Korabel and 
Barkai (201 1), in which the same question is answered using a random-walk model 
on a one-dimensional lattice. There, following the earlier work of Ovaskainen and 
Cornell (2003), they show that choosing how the random walk behaves at the inter- 
face of two regions of differing diffusion coefficient leads to different equilibrium 
probability densities for the system. Korabel and Barkai explain how to determine 
the correct interface behaviour of the random walk model using experimentally 
measurable quantities. Our approach differs in that, in our model, the interface 
behaviour is determined indirectly through the microscopic dynamics that we de- 
scribe. 

The results in Section [2] show that choosing a diffusion coefficient D(x) is not 
enough to fully specify diffusive dynamics in the absence of other assumptions. 
In particular, our results show that fixing D(x) is not enough to specify p eq (x), the 
equilibrium probability density for the position of the particle. In Section [3] we 
propose a framework for modelling state-dependent diffusion that makes this fact 
explicit. Rather than simply specifying a state-dependent diffusion coefficient, we 
specify diffusion coefficient D(x) and an equilibrium density p eq (x) which together 
with a detailed balance assumption (no-flux in equilibrium) completely determine 
the dynamics. We then introduce a new method for the numerical simulation of 
diffusive dynamics which makes use of our framework: the numerical method is 
expressed in terms of D(x) and p eq (x) and makes no reference to a drift term. The 
method consists of Euler-Maruyama steps for a purely diffusive Ito stochastic dif- 
ferential equation together with Metropolis rejections. The method is similar to the 
Metropolis-adjusted Langevin algorithm (Roberts and Tweedie, 1996; Bou-Rabee 
and Vanden Eijnden, 2010), except that there is no drift term in the Euler-Maruyama 
step and it is the Metropolis rejections that induce any drift in the trajectories. The 
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advantages of our method over the Euler-Maruyama method are that it samples 
space with the correct equilibrium density p eq (x) and performs well even when 
D(x) and p eq (x) have discontinuities. 

One area in which the proposed framework in Section [3] could be used is cellular 
biology. Although earlier models of chemistry in the cytoplasm of cells assumed 
that chemical species were 'well-mixed' and thus ignored diffusion, more recent 
models have taken the geometry of the cell and the diffusion coefficient of various 
molecules into account (Turner, Schnell, and Burrage, 2004). Effective diffusion 
coefficients of a molecule in a cell differ from that in water due to the crowding 
effect of other molecules. One approach is to model the motion of a molecule as 
diffusion with constant coefficient, but then have the effective diffusion be modified 
by interaction with other particles which are also included in the model (Ridgway, 
et al., 2008). Another is to not include the crowding particles in the model, but to 
model their effect with a modified diffusion coefficient (Hall and Hoshino, 2010). 
Given the inhomogeneity of the cytoplasm, we can expect that this effective diffu- 
sion coefficient of the molecule (and its equilibrium probability density) will vary 
with location within the cell. Our framework and numerical method are developed 
with this latter situation in mind. 

In Section [4] we conclude by explaining the relation between our results and three 
different interpretations of state-dependent diffusion: Ito, Stratonovich, and the 
Isothermal convention of Lancon, Batrouni, Lobry, and Ostrowski (2001). 



2 A Model System for State-Dependent Diffusion 



In the Introduction, we described a two-dimensional system consisting of a single 
particle diffusing inside a rectangular box and reflecting off the boundaries. The 
diffusion coefficient is twice as large on the right side of the box as the left. In this 
section, we demonstrate how to construct a family of deterministic Hamiltonian 
systems that approximates this behaviour on a coarse scale. 



In Subsection (2.1), we describe the random Lorentz gas (Dettmann, 2000), a deter- 
ministic system that when given a random initial condition yields constant-coefficient 



diffusion at a coarse scale. In Subsection (2.2), we show how to approximate the 
model system of the introduction by creating two adjacent domains of the random 



Lorentz gas within a bounding box. In Subsection (2.3), we describe numerical 
experiments with the box system demonstrating that the fraction of time a parti- 
cle spends on each side of the box cannot be determined solely from the values 



of the diffusion coefficient. In Subsection (2.4), we explain how the result of the 



numerical experiments in Subsection (2.3 ) can be predicted from properties of the 
dynamics of the microscopic system. 
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Fig. 2. Portions of the trajectory of the random Lorentz gas with = 0.5 at two levels of 
magnification, showing the discs and the trajectory of the moving particle. 

2. 1 Random Lorentz Gas 

Consider infinitely many discs with positions fixed in IR 2 . The centres of the discs 
are distributed randomly with uniform density subject to the constraint that the discs 
do not overlap. We consider a single point particle interacting with the discs in the 
following manner. Given an initial velocity and an initial position not on a disc, 
the particle moves with constant velocity until it meets a disc. Then it undergoes an 
instantaneous elastic reflection with the boundary of the disc, with the angle of inci- 
dence equalling the angle of reflection. This model is the two-dimensional random 
Lorentz gas (Dettmann, 2000), a mathematical formulation of a model originally 
due to Lorentz (Lorentz, 1905). We will always consider the case when the parti- 
cle has initial velocity of magnitude 1 . Figure [2] shows trajectories of the random 
Lorentz gas at two different scales. 

The random Lorentz gas has two parameters: the radius of the discs r, and the num- 
ber of discs per unit area A. The radius r can take any positive value. The density A 
has a maximum value =1/ (r 2 vT2) corresponding to the close-packed hexag- 
onal pattern of discs. We define to be the free volume fraction, the proportion of 
the area not occupied by discs. We have that = 1 — nr 2 X. The free volume can 
take any value in [0 m i n , 1) where m i n = 1 — n/\/V2 « 0.093, regardless of the value 
of r. The pair (r, 0) provides an alternative parametrization of the random Lorentz 
gas, with the advantage that is dimensionless. 

For = min adjacent discs are touching and the particle remains in a small region 
of the plane for its entire trajectory. For G (0 m i n , 1) the probability that two discs 
are touching anywhere in the plane is zero (Dettmann 2000, p. 325) and motion 
of the particle is conjectured to be diffusive (Dettmann and Cohen, 2000). (This 
contrasts with the situation of overlapping discs for which there is a phase transition 
to anomalous transport for dense enough placement of discs (Hofling, Munk, Frey, 



6 



and Franosch, 2008)). Specifically, if we start the particle at initial position x(0) not 
coincident with a disc and give it initial speed 1 and uniformly distributed direction 
in [0, lit) , then x(t), the position of the particle at time t is approximately distributed 
as a Gaussian random vector with mean and variance matrix 2DtI. Here / is the 
2x2 identity matrix and D denotes the diffusion coefficient. This conjecture is 
supported by analytical calculations (Dettmann, 2000; Ernst and Weyland, 1971) 
and numerical simulations (Dettmann and Cohen, 2000; Bruin, 1972). 

A stronger and more formal statement of the conjecture is stated in the language of 
weak convergence (Billingsley, 1999). Specifically, let x(t) denote the position of 
the moving particle at time t. Fix jc(0) to be some point not coincident with a disc 
and let the initial velocity be chosen as above. It is conjectured that 

x{nt)/y/n=> y/lD^Bit) 

as n — > oo, where B(t) is standard two-dimensional Brownian motion (meaning 
(B(t)) = 0, (B(t)B(t) T ) = I), D r A is the diffusion coefficient which depends on 
r and 0, and =>- denotes weak convergence in the space of continuous functions 
(Billingsley, 1999). Such a result holds for the certain periodic Lorentz gasses 
(Bunimovich and Sinai, 1980/81; Klages and Dellago, 2000), but remains open 
for the random Lorentz gas. (We chose for our study not to use the standard peri- 
odic Lorentz gas with discs centred on a hexagonal lattice since for large enough 
the particle undergoes superdiffusive motion in this case.) 

A scaling argument shows that, for fixed 0, D r A is proportional to r. To see this, 
letting | ■ | denote the Euclidean norm, observe that the coefficient D can be obtained 
as lim^ 00 (|x(f)| 2 )/4?, assuming x(0) = 0. If we rescale space by a factor R, we 
increase both the size of the discs and the distance the particle travels by a factor 
of R, without changing 0. So (|jc(?) 2 |) increases by a factor of R 2 .To maintain the 
speed of the particle as 1, we also have to rescale time, increasing t by a factor R. 
The net effect on the ratio (\x(t) \ 2 )/4t is to increase it by a factor R (Sanders 2005, 
Sec. 3.1.6). So 

D r ^=rf(<j>) (2.1) 

for some function / of . 

Figure [3] shows the relation between f(<j>) and that was computed using the tech- 
niques similar to those described in (Dettmann and Cohen, 2000). It appears that the 
function /(0) is continuous on its domain, it is monotonically increasing, f((j)) — > 
as — > 0nu n and /(0) goes to infinity as — > 1. Indeed, calculations from kinetic 
theory show that /(0) ~ 37r/[16(l — 0)] in the — > 1 limit (vanLeeuwen and Wei- 
jland, 1967; Bruin, 1972). 

Given any fixed diffusion coefficient D > there is a one-parameter family of 
choices ofr,0 such that D = D r ^ : for any e (0^,1), just choose r = 
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Fig. 3. The ratio of diffusion coefficient versus disc radius D/r as a function of free volume 
fraction for the random Lorentz gas. 

2.2 Box with Two Domains 

In order to investigate the main question of this paper, we take a rectangular box 
and divide it into two equal regions. Each region is filled with randomly placed 
discs, with different r and on each side. As in the Lorentz gas, the centres of 
the discs are placed uniformly at random with the condition that they not overlap 
with each other. Discs can intersect with the walls of the box but not the dividing 
line between the two sides of the box. Figures [4] and [5] show two examples of such 
configurations of discs. The dynamics of the point particle are the same as in the 
infinite random Lorentz gas, with the added condition that the particle reflects off 
the walls of the box. 

Suppose we are given diffusion coefficients D\,D2 > 0. We can choose ri,0i and 
r2,02 such that D\ = D ru ^ and £>2 = D ri ^ 2 . For small enough r\,r2 the dynamics 
of the particle in this system will be well- approximated by a particle that diffuses 
with coefficient D\ on the left side of the box and diffuses with coefficient D2 on the 
right side of the box. With appropriate choices of the parameters, we can investigate 
the question of the proportion of time the particle spends on each side of the box. 

2.3 Numerical Experiments 

We start the particle off at some position in the box not on a disc with velocity of 
magnitude 1 and randomly chosen direction. The motion of the particle is simulated 
with an event-driven simulation, computing a trajectory that is accurate up to the 
errors of floating point arithmetic (Dettmann and Cohen, 2000). Periodically the 
position of the particle is recorded. At the end of a long trajectory we compute the 
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Set-up 


n 


fa 


D X 




fa 


D 2 


Time on Right / Time on Left 


1 


0.3 


0.5 


0.09 


0.6 


0.5 


0.18 


0.99 


2 


0.09 


0.60 


0.047 


0.75 


0.30 


0.093 


0.50 



Table 1 

The parameters used in each side of the box in the two set-ups, and the ratio between 
amount of time spent by the particle on the right side of the box and the left side of the box 
in each set-up. 




Fig. 4. Set-up 1. Free volume fraction is the same on each side: <j>i = fa. Scatterer radius 
on the right is twice that on the left: 2r\ = r 2 , leading to 2D\ = D 2 . 

number of times the particle is on the right side of the box divided by the number 
of times the particle is on the left side of the box. 

We consider two choices of parameters on each side of the box, or set-ups, each 
of which leads to the effective diffusion coefficients satisfying Z>2 = 2D\. We have 
contrived the set-ups so that in Set-up 1, the Statistical Mechanics Prediction is 
correct, and that in Set-up 2, the Time-Change Prediction is correct. The parameters 
for each set-up are summarized in Table [T] 

In the first set-up 0i = fa = 0.5 and 2r\ = = 0.6. We show the position of the 
discs in the box in Figure |4j Over a trajectory of length 5 x 10 5 time units the ratio 
between the occupation times is approximately 1, as we show in the table above. 
This result agrees with the Statistical Mechanics Prediction. 

In the second set-up 202 = fa = 0.60 and 8.3ri ss r2 = 0.75. We show the position 
of the discs in the box in Figure [5} Over a trajectory of length 7 x 10 6 time units the 
ratio between the occupation times is approximately 1/2, as we show in the table 
above. This result agrees with the Time-Change Prediction. 

Thus, by fixing the parameters appropriately, both the Statistical Mechanics Pre- 
diction and the Time-Change Prediction can be seen to be correct for the given 
mesoscopic behaviour. 
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Fig. 5. Set-up 2. Free volume fraction is twice as big on the left: 0i = 202- Scatterer radii 
are chosen so that 2D\ = D2. 

2.4 Analysis and Discussion 

We first explain that, given the properties of random Lorentz gas, the results of the 
simulation in the previous subsection are predictable. We first note that the system 
is ergodic. The ergodicity of periodic Lorentz gas is shown by Sinai (1970), since 
it is equivalent to dispersing billiards on the torus. Our system is not dispersing 
because the walls of the box are not convex, but Sinai (1970, Sec. 9) shows how 
ergodicity still holds in this case by unfolding the box to obtain a periodic domain. 

Ergodicity implies that the duration of time that the particle spends in a region of 
phase space is proportional to the volume of the region. For our system, this implies 
that the amount of time spent by the particle on each side is proportional to the free 
volume fraction on each side. For fixed 0, the parameters r and D are irrelevant to 
the proportion of time the particle spends on each side of the box. In Set-up 1 above, 
01 = 02 so the ratio of times spent on each side of the box is equal and the Statistical 
Mechanics Prediction is correct. In Set-up 2, 0i = 202, and so the particle spends 
twice as much time on the left side of the box, and the Time-Change Prediction is 
correct. 

Despite appearances, the principles of statistical mechanics are not violated in Set- 
up 2. There are two ways to reconcile the apparent contradiction. Firstly, one can 
say that at the microscopic scale statistical mechanics is not violated because there 
are not an equal number of states on each side. The number of states is propor- 
tional to the free volume fraction on both sides, and so the system spends more 
time where the free volume fraction is higher. The other way of reconciling the 
disagreement is at the mesoscale. Suppose we divide the box into many rectangular 
cells of equal area, each much smaller than the whole box, but much larger than 
the size of the discs. Each cell corresponds to a mesoscopic state which the parti- 
cle may be in. Since the system is in the microcanonical ensemble (no exchange 
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of energy with the outside) the probability of the system in equilibrium being in 
a particular mesoscopic state is determined by the entropy of that state, where the 
entropy of a mesoscopic state is proportional to the logarithm of the amount of mi- 
croscopic states it contains. The greater free volume fraction on the left in Set-up 
2, implies that mesoscopic states have greater entropy there, and so the system will 
spend more time on the left side of the box. 

On the other hand, in Set-up 1, the Time-Change Prediction proves to be wrong. 
This means that the motion of the particle in the box is not well-described by the 
drift- free Ito stochastic differential equation ( |1.1| ). Since, by the properties of the 
uniform random Lorentz gas, < | 1 . 1[ ) is a good model for the dynamics of the particle 
within each of the two regions of constant disc radius r, it must be that the equation 
is no longer a good model at the boundary of the two regions. 

We see that there are choices of 0i,n and 02,^2 sucn that either the Time-Change 
Prediction or the Statistical Mechanics Prediction are correct, while still D2 = 2D\ . 
Indeed, for any D\ and D2, parameters can be chosen to induce arbitrary ratios 
between the times spent on the left-hand side and the right-hand side. 

Although for generic values of the parameters in our Lorentz gas model neither pre- 
diction will be valid, we point out that the Statistical Mechanics Prediction holds 
for a natural set of parameter settings, whereas the same is not true of the Time- 
Change Prediction. For the Time-Change Prediction to be correct, it is necessary 
for the ratio between <j>\ and 02 to match the ratio between D2 and D\. We see 
no natural way that the parameters in our model may be set for this matching to 
occur. On the other hand, for the Statistical Mechanics Prediction to be correct it 
is necessary that 0i and be equal. There is at least one case where this condi- 
tion approximately holds for a naturally occurring system. Consider a situation in 
which both 0i,02 ~ 1- This requires no fine tuning, only that the discs take up a 
small fraction of the total volume. Choosing r\ and to be unequal leads to differ- 
ent diffusion coefficients on each side, but the particle still spends approximately 
equal proportions of time in each region. Likewise, in any physical system where 
a particle diffuses by interacting with small, sparsely placed scatterers, we expect 
the Statistical Mechanics Prediction to be correct. 



3 A Proposal for modelling with state-dependent diffusion 

The numerical experiments of the previous section demonstrate that the equilibrium 
density p eq (x) is not determined solely by the local diffusion rate D(x), even in 
situations with no external forces acting on the particle. This raises the practical 
issue of how to model systems with state-dependent diffusion. Ideally, mesoscopic 
diffusive models would be derived from microscopic models via an asymptotic 
technique such as the van Kampen system-size expansion (van Kampen 2007, Ch. 
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XI.3). But in many circumstances, deriving a realistic microscopic model may be 
impractical. 



Instead, in Subsection ( |3.1[ ) we describe a more phenomenological approach. We 
assume that the modeller posits an isotropic state-dependent diffusion rate D(x) 
and an equilibrium density p eq (x). We then derive a drift coefficient a(x) for an Ito 
stochastic differential equation with diffusion coefficient D(x) that gives the desired 
p e q{x). In Subsection (3.2) we show how the Lorentz gas system of Section[2]can 
be modeled at a mesoscopic level in this way. In Subsection (3.3 ) we then describe 
an algorithm for simulating such stochastic differential equations (SDEs) that does 
not refer to a(x) but is expressed in terms only of D(x) and p eq (x). 



3.1 A modelling framework for state-dependent diffusion 



We model the diffusion in k spatial dimensions by the Ito SDE 



dX(t) = a(X)dt + y/2D(X)dB(t), (3.1) 

where we will determine a(x) in terms of p eq (x) and D(x). Here B[t) is standard 
^-dimensional Brownian motion. (In alternate notation, we write this equation as 
dX(t)/dt = a(x) + a/2D(X)t7 (t) where Tj is fc-dimensional Gaussian white noise.) 
The Fokker-Planck equation for this system is 

J^pfoO = - V- [a(x)p{x,t)} +A[D(x)p(x,t)} = -V-J(x) 

where J is the probability flux. Since the problem is underdetermined as stated, 
we will stipulate that the probability flux vanishes in equilibrium, which is the 
same as detailed balance holding; see van Kampen (2007, Ch. XI.4) for criteria on 
systems under which this condition holds. We choose a(x) so that J(x) is zero in 
equilibrium: 

J(x) = a(x)p eq (x) - V[D(x)peq(x)] = 0. 
Solving for a(x) gives 

a(x) = — — V(D(jc)peq(jc)) = VD(x)+D(x)V]np eq (x). (3.2) 

Peq (X) 

Thus given an equilibrium density p eq (x) and diffusion coefficient D(x) the appro- 
priate Ito SDE is 

dX(t) = (VD(X)+D(X)V\np eq (X))dt+ v / 2D(X)dB(t). (3.3) 



The Fokker-Planck equation of ( |3.3[ ) is 

-p(x,t)=V-[-V(D(x)p eq (x))(p(x)/p sq (x)) + V(D(x)p(x))} 
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or simplifying, 

j t PM = V-[D(x)p eq (x)V(p(x)/p eq (x))}, (3.4) 

which is the special isotropic case of Equation XI.4.14 of van Kampen (2007). 
Note that if p eq (x) is constant with respect to x, which corresponds to the Statistical 
Mechanics Prediction being true, then a(x) = VD(x) and ( |3.3| ) reduces to 

dX(t) = VD(x)dt+y/2D(x)dB(t). 

On the other hand, to obtain a drift-free Ito SDE ((a(x) = 0) in this framework 
requires D(x)p cq (x) to be constant in x, which means that p eq (x) is determined 
completely by D(x). 



3.2 Connection to the Lorentz gas model 



We explain the connection between the framework of Subsection (3.1 ) and the mi- 
croscopic Lorentz gas model of Section[2} Given an instance of our random Lorentz 
gas model with two domains, what are the corresponding functions D(x), p eq (x) in 



( |3.4[ ), the mesoscopic equation for p? 



Suppose on the left side of the box the Lorentz gas model has parameters r\ and 
01 and on the right it has parameters r2 and 02- The diffusion coefficients on each 



side, D{ and Z>2 are determined by the relation ( |2.1[ ). To determine the equilibrium 
probability density of the particle on each side, let 2A be the total area of the box, 
so that each side has area A. We know that p eq ; , the probability density on side i, is 
proportional to 0,-, and thus p eq .i/Peq,2 = 0i/02- We also know that since the total 
probability must be 1, p eq ,\A + p eq .2^ = 1. Solving for p eq j gives 



y eq,( 



A(0i+0 2 ) 



for i = 1,2. We define D(x) and p eq (x) for x = (x\ 1 X2) in the box by 



D(x) 



Di, for x\ < 0, 
Z>2, for x\ > 0, 



PeqO) 



Peq,!, forXl < 0, 

p eqi 2, for^i > 0. 



These functions determine a(x), the drift coefficient in ( |3.1| ), via ( 3/1) . The drift 
a(x) is zero everywhere except along the line x\ = where it is not defined. 

We determine the appropriate boundary conditions for p along the boundary line 



x\ = 0. In order for the right hand side of (3.4) to be well defined we require 
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p(x) I Peq(jc) to be continuous. So for any x along the boundary line we must have 

P(*~) = P(* + ) 
pe q (*~) p eq (x+) 

where x~ denotes taking the limit from the left, and x + denotes taking the limit 
from the right. For our particular choice of p eq this gives 

P(*~) _ P(* + ) 

Peq,l Peq,2 

The second boundary condition comes from assuming continuous flux across the 
boundary line: 

D(X p^ X W— H \ \ =D(x + )p eq (x + ) 1 r- H \ ' 

J dx x p eq (x-) ^ qv J dx l p eq (x+) 

For our particular choices of D and p eq , since p eq is constant away from the line 
x\ = 0, this gives the boundary conditions 

A possible direction for further investigation is to consider Lorentz gas models 
where disc radius r and free volume fraction vary smoothly with x, and to deter- 
mine what D(x) and p eq (*), an d hence a(x) are in this case. 



3.3 Numerical Simulation 



If both p eq and D are smooth then the Euler-Maruyama scheme or Milstein's method 



is effective for simulating ( |3.3[ ) (see Higham 2001, for example). If either D(x) or 
p eq (x) are discontinuous with respect to x, then the drift term will have a singular- 
ity that is not resolved by standard time-integration schemes, and numerical sim- 
ulations may not correctly approximate the equilibrium density. We propose using 
a variant of the Metropolis-adjusted Langevin algorithm (MALA) as introduced 
by Roberts and Tweedie (1996) and analyzed by Bou-Rabee and Vanden-Eijnden 
(2010). MALA is obtained by taking a convergent numerical method for the SDE 
(in this case, Euler-Maruyama) and introducing Metropolis step-rejections in order 
to have a discrete-time process with the correct equilibrium density. Our approach 
is different in that we start with a convergent method for only the diffusive part of 
the Ito SDE, and then we introduce step rejections to induce the correct drift and 
equilibrium density. 

Let h be the step length of our numerical discretisation and let X n be the numerical 



approximation to X(nh), where X is the solution to the Ito SDE ( |3.3| ) for n > 0. 
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Given a numerical value X n , we let the trial value for the next step be defined by 

X* +1 =X n + y/2D(X n )[B((n + l)h)- B(nh)} (3.5) 
where B is standard ^-dimensional Brownian motion. Then X n+ \ is given by 

X n +i = \ (3.6) 
X„ otherwise, 



where 



a h (x,y) = mm 1, — — — ^— . 



and %k,k > 1 is an independent, identically distributed sequence of random vari- 
ables, uniform on [0, 1] and independent of B. Here 



qh(x ' y) = 7wm e 



-(x-y) 2 /4hD(x) 



is the transition probability density for X* +l being at y given that X n is at x. The 



definition of X* +l in p.5[ ) is just the Euler-Maruyama method for the Ito stochastic 



differential equation dX = ^2D(X)dB. The Metropolis rejection procedure ( |3.6| ) 



accepts the Euler-Maruyama step with probability oc^(Z n ,X* +1 ) and rejects it oth- 
erwise. 

The Metropolis rejection procedure guarantees that the process X n ,n> 1 has p eq (x) 
as its equilibrium density. On the other hand, in any region of the state space where 
D(x) and p eq (x) are constant with respect to x, ah is 1 and so the method reduces 
to the Euler-Maruyama method for the constant coefficient diffusion without drift. 
Future work will investigate rigorously the convergence of the above method to the 



solutions of (3.3). 



Here we numerically demonstrate the convergence of the method for the simple 
case where 

D(x) = 

and 

' l if* e [-1,1] 

-eq**) = \ 

otherwise. 




Equation ( |3.3[ ) with this choice of D(x) and p sq (x) models a particle diffusing on the 
interval [—1,1] with reflecting boundary conditions at ± 1 and a piecewise constant 
diffusion coefficient. The reflecting boundary conditions are conveniently imple- 
mented by our choice of p eq (x). 
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-8— t— 9--B- 



Fig. 6. Equilibrium density p eq and local diffusion coefficient D for method ( |3.5| > applied to 
a simple one-dimensional SDE. Shown are results for h = 0.01 (dotted), h = 0.001 (dashed), 
and h = 0.0001 (solid), along with the exact values for the SDE. 



Figureoshows the results of simulating ( |3.3| ) with these choices of D(x) and p eq (x 



using the method we have described. To show results, we divide the interval [—1,1] 
into 20 equal subintervals, and plot the density for the amount of time the particle 
spends in each subinterval. We also plot the effective diffusion coefficient for each 
bin, which we define to the the average observed value of (X n+ \ —X n ) 2 /2h over the 
trajectory, for all n such that X n lies in the given bin. Simulations were conducted 
with h = 0.01 , 0.001 , 0.0001 and for trajectories long enough so that standard statis- 
tical errors in the plots are smaller than the symbols used. We see that for all values 
of h the equilibrium density p eq (*) is correctly reproduced. The effective diffusion 
coefficient converges to D as h goes to zero. 



4 Discussion: Ito, Stratonovich, or Isothermal 



We conclude by discussing our results in the context of the apparent ambiguity be- 
tween Ito, Stratonovich, and Isothermal interpretations of stochastic integrals (Lau 
and Lubensky, 2007; Volpe, et al., 2010; Kupferman, et al., 2004). For the purposes 
of discussion, we consider a particle moving in one spatial dimension whose po- 
sition at time t is X(t). We assume that X(t) is a Markov stochastic process with 
continuous sample paths. We model the motion of the particle with the stochastic 
differential equation (SDE) 

dX(t)=a(X(t))dt + b(X(t))dB (4.1) 

where B(t) is standard Brownian motion. We call a(x) the drift and b(x) the diffu- 
sion of the SDE. 

As is well known (van Kampen, 2007; Gardiner, 2004), unless we specify a par- 



ticular interpretation, the SDE (4.1 ) does not unambiguously define the stochastic 
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process X(t). To see this, we integrate ( |4~T 



I over [0, T] to get 



X(t)-X(0)= a(X(t))dt+ b(X(t))dB(t). 
Jo Jo 

The first term on the right has a unique interpretation as a Riemann integral, but the 
second term cannot be simply viewed as a Riemann-Stieljes integral, since B is not 
of bounded variation. If we compute the second term as the limit of Riemann sums, 
the answer depends on where in each subinterval the argument b(X(t)) is evaluated. 
For example, choosing h = T/N, t n = hn and letting B n = B[t n ) and X n = X(t n ), 
suppose we take the integral with respect to B to be 

f T b{X(t))dB(t) = \mxfb{X* n ){B n+x -B H ), 
Jo ^°„= 

where 

X* = (l-a)X n + aX n+l . 

Famously, unless b(x) is a constant, the limit depends on the choice of a (Volpe, 
et al., 2010). If we choose a = 0, we obtain the Ito interpretation of the integral, 
which yields a stochastic process X(t) with Fokker-Planck equation 

j f p(x,t) = - ~ [a(x)p(x,t)] + \^2 [K*) 2 P(*,0] ■ 

If we choose a = 1/2, we obtain the Stratonovich interpretation of the integral, 
which yields a process X(t) with Fokker-Planck equation 



d . . d 
dt pM = -dx 



a(x)p(x,t) + -b(x)b'(x)p(x,t) 



+ 



2dx 2 



[b{x) 2 p{x,t 



Note that the Fokker-Planck equation shows that the Stratonovich interpretation 
of the SDE with drift a(x) and diffusion b(x) yields the same stochastic process 
as the Ito interpretation of the SDE with drift a(x) +b(x)b'(x)/2 and diffusion 
b(x) (Gardiner 2004, p. 99). Finally, if we choose a = 1 we obtain the anti-Ito 
or Isothermal interpretation (Lau and Lubensky, 2007; Volpe, et al., 2010), which 
yields a process with Fokker-Planck equation 



j^pfoO = [a{x)p(x,t)+b{x)b'(x)p(x,t)] [b(x) 2 p(x,t)] . 

In this case, the Isothermal interpretation of the SDE with drift a(x) and diffusion 
b(x) gives the same stochastic process as the Ito interpretation of the SDE with drift 
a(x) +b(x)b'(x) and diffusion b(x) (Lau and Lubensky, 2007). 

As we can see in the various Fokker-Planck equations above, if we fix a(x) and 
b(x), varying the parameter a gives different stochastic processes for the motion of 
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a 


Interpretation of 
oiocndsiic integral 


Fokker-Plank Equa- 
tion tor a = U 


Special Properties 
when a = 





Ito 




(X(f)--X(O)) =0 for all f. 


1/2 


Stratonovich 


&> = *&{*&[*P]} 




1 


Isothermal 


|p = ^{^p} 


p eq = COnSt. 



Table 2 

Summary of some properties of the Ito, Stratonovich, and Isothermal interpretations of 



the particle. This fact may make it seem like there should be a physically correct 
choice of the parameter a. We argue that this is false. For any fixed a the range 
of stochastic processes that can be captured by an appropriate choice of a(x) and 
b(x) is the same. For example, suppose we fix a choice of a(x) and b(x) and choose 
to interpret ( |4.l[ ) with a given a £ [0, l]. The process defined is identical to what 
we would obtain with the Ito interpretation (a = 0) of the SDE with drift a(x) + 
ab'(x)b(x) and diffusion b(x). 



Though the families of stochastic processes described using each convention are 
the same, it may still be the case that some choice of a is more natural or conve- 
nient for some purposes than others. Frequently, the rationale is based on the idea 
that if the drift is zero, then X(t) should have certain properties. For example, if one 
wants (X(t) — X(0)) = for all t when a(x) = 0, regardless of b(x), then the Ito 
convention with a = guarantees this. If one wants that when a(x) = the equi- 
librium density is constant, then the Isothermal convention with a = l guarantees 



this. We summarize the properties of the various interpretations of the SDE ( |4.l| ) 
when a(x) =0 in Table [2j 

The question we posed in the Introduction may be rephrased as follows: in the 
absence of external forces, and given a diffusion b(x) = v/2Z)(x), what is the cor- 
rect choice of parameter a and drift a(x) to model the motion of the particle? A 
natural way to approach the problem is to interpret the absence of external forces 
as meaning that a(x) = 0. Then the problem boils down to the choice of a: the 
Statistical-Mechanics Prediction follows from taking a = l, and the Time-Change 
Prediction follows from taking a = 0. The results in Section [2] showed that neither 
answer is justified universally. 

In Section [3] we recommended a different approach. We fix a and then choose a(x) 
to generate the desired equilibrium distribution. As we have explained here, the 
choice of a is not crucial once we allow a non-zero a(x). Accordingly, we have 
chosen a = 0, corresponding to Ito calculus. This is the main choice in the mathe- 
matics literature, and numerical methods such as the Euler-Maruyama method take 
a particularly simple form with it. Once we have made this choice of a, we are 
free to choose a(x) appropriately. In Section [3] we chose a(x) to ensure a given 
equilibrium density. 
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Beyond the particular needs of the present work, we believe the framework we 
describe in Section [3] provides a natural and flexible way to model diffusive sys- 
tems. In situations where a researcher is confident for physical reasons that the 
equilibrium probability is constant, then our framework takes a simple form. The 
numerical method we present functions even when the diffusion and equilibrium 
density are discontinuous. Future work will study the convergence properties of the 
method, as well exploring applications of our general framework. 
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